
/*
REPLICATION OF FIGURES FOR PAPER: "THE RISE AND PERSISTENCE OF ILLEGAL CROPS: EVIDENCE FROM A NAIVE POLICY ANNOUNCEMENT"
Authors: Mounu Prem, Juan F. Vargas, and Daniel Mejía
*/

cls	
clear all 
set more off
set seed 11110


** Path to main folder
	global path "~/Dropbox/Research/Research Projects/Joint with Juan/Coca/Paper Replication/"
	
** Subfolders
	global Data 	"${path}Data/"
	global Figures 	"${path}Figures/"
	global Tables 	"${path}Tables/"

** Set figures folder
	cd "${Figures}"	

*****************************************************
** MAIN PAPER FIGURES
*****************************************************	

************
// Figure 1. Evolution of coca cultivation
************
	* Load dataset
	use "${Data}Data_CountryCoca.dta", clear
	
	two (connect colombia year , msize(small) msymbol(square) mcolor(black) lcolor(black)) ///
		(connect peru year , msize(small) msymbol(square) mcolor(gs5) lcolor(gs5)) ///
		(connect bolivia year , msize(small) msymbol(square) mcolor(gs10) lcolor(gs10)) , scale(1.2) xlabel(1999(3)2018) ///
		ytitle("Coca cultivation per 1,000 hectares") legend(label(1 "Colombia") label(2 "Bolivia") label(3 "Peru"))  graphregion(fcolor(white)  /// 
		lcolor(white) ifcolor(white) ilcolor(white)) plotregion(fcolor(white) lcolor(white) ifcolor(white) ilcolor(white))
		
	gr export "Figure_1.pdf", replace
	
************
// Figure 2. Raw data dynamics and point estimates
************

	** Panel A. Raw data: Suitability
	* Load dataset
	use "${Data}MainData.dta", clear

	local treat "suit_median"
	gen var_Coca	= coca_area_grid if `treat'==1 
	gen var_nonCoca	= coca_area_grid if `treat'==0
								
	collapse (mean) var_Coca var_nonCoca , by(year) 
	
	replace var_Coca=. if var_Coca==0
	replace var_nonCoca=. if var_nonCoca==0
	
	tsset year, y
	
	twoway  (tsline var_Coca, lcolor(black) lpattern(dash)  lwidt(thick)) ///
		(tsline var_nonCoca, lcolor(gs5) lpattern(solid)  lwidt(thick)), ///  
		xline(2013, lp(dash_dot) lcolor(gs10)) ///
		xline(2016.9, lp(dash_dot) lcolor(gs10)) ///
		graphregion(fcolor(white) lcolor(white) ifcolor(white) ilcolor(white)) ///
		xtitle(" ") xlabel(2011(1)2018, angle(ninety) labsize(vsmall)) ///					
		legend(label(1 "High Suitability") label(2 "Low Suitability"))  
		*legend(label(1 "High probability") label(2 "Low probability")) 	
	
	gr export "Figure_2A.pdf", replace
	
	** Panel B. Raw data: Probability of substitution
	* Load dataset
	use "${Data}MainData.dta", clear

	local treat "pr_PNIS_50"
	gen var_Coca	= coca_area_grid if `treat'==1 
	gen var_nonCoca	= coca_area_grid if `treat'==0
								
	collapse (mean) var_Coca var_nonCoca , by(year) 
	
	replace var_Coca=. if var_Coca==0
	replace var_nonCoca=. if var_nonCoca==0
	
	tsset year, y
	
	twoway  (tsline var_Coca, lcolor(black) lpattern(dash)  lwidt(thick)) ///
		(tsline var_nonCoca, lcolor(gs5) lpattern(solid)  lwidt(thick)), ///  
		xline(2013, lp(dash_dot) lcolor(gs10)) ///
		xline(2016.9, lp(dash_dot) lcolor(gs10)) ///
		graphregion(fcolor(white) lcolor(white) ifcolor(white) ilcolor(white)) ///
		xtitle(" ") xlabel(2011(1)2018, angle(ninety) labsize(vsmall)) ///					
		legend(label(1 "High probability") label(2 "Low probability")) 	
	
	gr export "Figure_2B.pdf", replace


	** Panel C. Point estimates: Suitability
	* Load dataset
	use "${Data}MainData.dta", clear

	foreach i of num 2011(1)2018 {		
	    gen y`i' = (year==`i')		
    }
	
	foreach i of num 2011(1)2018 {	
    	gen suit_y`i' = suitability * y`i'	
    }
	
	global I "suit_y2011 suit_y2012 suit_y2014 suit_y2015 suit_y2016 suit_y2017 suit_y2018"
	local parm_v = "suit_y2013"
	
	reghdfe coca_area_grid $I , absorb(muni_code i.coddepto#i.year) vce(cluster muni_code)	
		* Plot
		preserve
			estimate store coef 
			parmest, norestore
			rename t tstat
			g year = 2010 + _n
			replace year = year + 1 if _n > 2
			set obs 8
			replace year = 2013 if _n == _N
			replace estimate = 0 if year == 2013
			replace min95 = 0 if year == 2013
			replace max95 = 0 if year == 2013
			replace parm="`parm_v'" if parm==""
			sort year
			
			twoway (connected estimate year , lc(gs0) ) ///
			(rcap min95 max95 year , lc(gs12) fc(gs12)), ///
			xlabel(2011(1)2018) ///
			xtitle("Years") ///
			ytitle("Coefficient") ///
			yline(0, lc(red) ) ///
			ylabel(, grid ) ///
			xline(2016.9, lp(dash_dot) lcolor(ebblue)) ///
			graphregion(color(white)) /// 
			legend(off)
		restore
	
		gr export "Figure_2C.pdf", replace

		
		
	** Panel D. Point estimates: Probability of substitution
	* Load dataset
	use "${Data}MainData.dta", clear

	foreach i of num 2011(1)2018 {		
	    gen y`i' = (year==`i')		
    }
	
	foreach i of num 2011(1)2018 {	
    	gen pr_PNIS_y`i' = pr_PNIS * y`i'	
    }
	global I "pr_PNIS_y2011 pr_PNIS_y2012 pr_PNIS_y2014 pr_PNIS_y2015 pr_PNIS_y2016 pr_PNIS_y2017 pr_PNIS_y2018"
	local parm_v = "pr_PNIS_y2013"
	
	reghdfe coca_area_grid $I , absorb(muni_code i.coddepto#i.year) vce(cluster muni_code)	
		* Plot
		preserve
			estimate store coef 
			parmest, norestore
			rename t tstat
			g year = 2010 + _n
			replace year = year + 1 if _n > 2
			set obs 8
			replace year = 2013 if _n == _N
			replace estimate = 0 if year == 2013
			replace min95 = 0 if year == 2013
			replace max95 = 0 if year == 2013
			replace parm="`parm_v'" if parm==""
			sort year
			
			twoway (connected estimate year , lc(gs0) ) ///
			(rcap min95 max95 year , lc(gs12) fc(gs12)), ///
			xlabel(2011(1)2018) ///
			xtitle("Years") ///
			ytitle("Coefficient") ///
			yline(0, lc(red) ) ///
			ylabel(, grid ) ///
			xline(2016.9, lp(dash_dot) lcolor(ebblue)) ///
			graphregion(color(white)) /// 
			legend(off)
		restore
	
		gr export "Figure_2D.pdf", replace

	******* FIX ME	
	** Panel E. Larger sample: All grids
	/*
	global T "d_3 d_2 d0 d1 d2 d3 d4"
	
	reghdfe COCA_C_ $T, abs(OBJECTID year year#muni_code) vce(cl muni_code)
		* Plot
		preserve
			estimate store coef 
			parmest, norestore
			keep if _n <= 7
			rename t tstat
			
			g year = 2010 + _n
			replace year = year + 1 if _n > 2
			set obs 8
			replace year = 2013 if _n == _N
			replace estimate = 0 if year == 2013
			replace min95 = 0 if year == 2013
			replace max95 = 0 if year == 2013
			replace parm="`parm_v'" if parm==""
			sort year
			
			twoway (connected estimate year , lc(gs0) ) ///
			(rcap min95 max95 year , lc(gs12) fc(gs12)), ///
			xlabel(2011(1)2018) ///
			xtitle("Years") ///
			ytitle("Coefficient") ///
			yline(0, lc(red) ) ///
			ylabel(, grid ) ///
			graphregion(color(white)) /// 
			legend(off)
			gr export "fig_2_coca_grid_rf.pdf", replace
		restore
	*/	
		
	** Panel F. Larger sample: Within coca grids
	/*
	global T "d_3 d_2 d0 d1 d2 d3 d4"
	
	reghdfe COCA_C_ $T if ever_coca == 1, abs(OBJECTID year year#muni_code) vce(cl muni_code)
		* Plot
		preserve
			estimate store coef 
			parmest, norestore
			keep if _n <= 7
			rename t tstat
			
			g year = 2010 + _n
			replace year = year + 1 if _n > 2
			set obs 8
			replace year = 2013 if _n == _N
			replace estimate = 0 if year == 2013
			replace min95 = 0 if year == 2013
			replace max95 = 0 if year == 2013
			replace parm="`parm_v'" if parm==""
			sort year
			
			twoway (connected estimate year , lc(gs0) ) ///
			(rcap min95 max95 year , lc(gs12) fc(gs12)), ///
			xlabel(2011(1)2018) ///
			xtitle("Years") ///
			ytitle("Coefficient") ///
			yline(0, lc(red) ) ///
			ylabel(, grid ) ///
			graphregion(color(white)) /// 
			legend(off)
			gr export "fig_2_coca_grid_rf.pdf", replace
		restore
	*/			
	
************
// Figure 3. Change in coca production, coca suitability, and substitution probability
************
	* See Maps.R code

*****************************************************
** APPENDIX FIGURES
*****************************************************	

************
// Figure A.1. Main milestones and announcements during the peace negotiations
************
	* Hand made
	
	
************
// Figure A.2. Dynamic specification in a longer series 2008 to 2018
************
	* Load dataset
	use "${Data}MainData_Long.dta", clear
	
	sum t
	local min = -r(min)
	forval i = 1/`min' {
		g d_`i'_p = (t== -`i') * pr_PNIS
		g d_`i'_s = (t== -`i') * suitability		
	}
	
	sum t
	forval i = 0/`r(max)' {
		g d`i'_p = (t== `i') * pr_PNIS
		g d`i'_s = (t== `i') * suitability		
	}
	
	global Tp 	"d_6_p d_5_p d_4_p d_3_p d_2_p d0_p d1_p d2_p d3_p d4_p"
	global Ts 	"d_6_s d_5_s d_4_s d_3_s d_2_s d0_s d1_s d2_s d3_s d4_s"
	
	reghdfe sh_coca $Ts , absorb(muni_code year i.coddepto#i.year) vce(cluster muni_code) 
		preserve
			estimate store coef 
			parmest, norestore
			rename t tstat
			g year = 2007 + _n
			replace year = year + 1 if _n > 5
			set obs 11
			replace year = 2013 if _n == _N
			replace estimate = 0 if year == 2013
			replace min95 = 0 if year == 2013
			replace max95 = 0 if year == 2013
			sort year
			
			twoway (connected estimate year , lc(gs0) ) ///
			(rcap min95 max95 year , lc(gs12) fc(gs12)), ///
			xlabel(2008(2)2018) ///
			xtitle("Years") ///
			ytitle("Coefficient") ///
			yline(0, lc(red) ) ///
			ylabel(, grid ) ///
			graphregion(color(white)) /// 
			legend(off)
			gr export "AppendixFigure_2A.pdf", replace
		restore
	
	reghdfe sh_coca $Tp , absorb(muni_code year i.coddepto#i.year) vce(cluster muni_code) 
		preserve
			estimate store coef 
			parmest, norestore
			rename t tstat
			g year = 2007 + _n
			replace year = year + 1 if _n > 5
			set obs 11
			replace year = 2013 if _n == _N
			replace estimate = 0 if year == 2013
			replace min95 = 0 if year == 2013
			replace max95 = 0 if year == 2013
			sort year
			
			twoway (connected estimate year , lc(gs0) ) ///
			(rcap min95 max95 year , lc(gs12) fc(gs12)), ///
			xlabel(2008(2)2018) ///
			xtitle("Years") ///
			ytitle("Coefficient") ///
			yline(0, lc(red) ) ///
			ylabel(, grid ) ///
			graphregion(color(white)) /// 
			legend(off)
			gr export "AppendixFigure_2B.pdf", replace
		restore	
		
		

************
// Figure A.3. Permutation test
************
	* Load dataset
	use "${Data}MainData.dta", clear

	global B 500
	preserve
		clear 
		g n = .
		g b1 = .
		g b2 = .
		
		save "${Data}PlaceboInference.dta", replace
	restore

	egen id = group(muni_code)
	
	forval b = 1 /$B {	
	
		preserve
			keep suitability 
			duplicates drop
			g aux = runiform(0,1)
			sort aux
			rename suitability placebo1
			g id = _n
			tempfile temp1
			save `temp1'.dta, replace
		restore
		
		
		preserve
			keep pr_PNIS
			duplicates drop
			g aux = runiform(0,1)
			sort aux
			rename pr_PNIS placebo2
			g id = _n
			tempfile temp2
			save `temp2'.dta, replace
		restore
		
		
		merge n:1 id using `temp1'.dta
		drop _m

		merge n:1 id using `temp2'.dta
		drop _m

		g placebo_post1 = placebo1 * announcement
		g placebo_post2 = placebo2 * announcement
		
		* Main: Dummy
		reghdfe coca_area_grid placebo_post1 , absorb(muni_code year i.coddepto#i.year) vce(cluster muni_code)		
		g b1 = _b[placebo_post1]
		
		reghdfe coca_area_grid placebo_post2 , absorb(muni_code year i.coddepto#i.year) vce(cluster muni_code)		
		g b2 = _b[placebo_post2]
		g n = `b'
		
		preserve
			keep b1 b2 n 	
			duplicates drop
			
			append using "${Data}PlaceboInference.dta"
			
			compress
			save "${Data}PlaceboInference.dta", replace
		restore
		
		drop placebo_post1 placebo_post2 placebo2 b1 b2 n
	}
	
	use "${Data}PlaceboInference.dta", clear
	
	kdensity b1, ///
		xline(0.35, lc(gray)) ///
		graphregion(fcolor(white) lcolor(white) ifcolor(white) ilcolor(white)) ///
		plotregion(fcolor(white) lcolor(white) ifcolor(white) ilcolor(white)) scale(1.2) xlabel(-0.2(0.1)0.4) ///
		legend(off) title("") ///
		xtitle("Coefficient") 
	
	graph export "AppendixFigure_3A.pdf"	, replace

	
	kdensity b2, ///
		xline(26.48, lc(gray)) ///
		graphregion(fcolor(white) lcolor(white) ifcolor(white) ilcolor(white)) ///
		plotregion(fcolor(white) lcolor(white) ifcolor(white) ilcolor(white)) scale(1.2) xlabel(-2(4)28) ///
		legend(off) title("") ///
		xtitle("Coefficient") 	
	
	graph export "AppendixFigure_3B.pdf"	, replace
	
************
// Figure A.4. Change in conflict-related events by coca suitability and substitution probability
************

	* Load dataset
	use "${Data}MainData.dta", clear
	
	collapse viol_oagv2d_* att_oagv2d_* suitability pr_PNIS indrural discapital IPM lpobl_tot suit_median pr_PNIS_50 (firstnm) coddepto, by(muni_code)

	local i "viol_oagv2d"
	g `i' =  `i'_18 -`i'_13
	replace `i' = 0 if `i' == .
	
	** Panel A. Violence: Suitability
	label def xx2 0 "Low Suitability" 1 "High Suitability" 
	label values suit_median xx2
	graph bar `i', over(suit_median) ytitle("") ///
		graphregion(fcolor(white))  bar(1, fcolor(black)) bar(2, fcolor(black))
		
	gr export "AppendixFigure_4A.pdf", replace
	
	** Panel B. Violence: Probability of substitution
	label var pr_PNIS_50 "Probability of substitution"
	label def xx 0 "Low probability" 1 "High probability"
	label values pr_PNIS_50 xx
	graph bar `i', over(pr_PNIS_50) ytitle("") ///
		graphregion(fcolor(white))  bar(1, fcolor(black)) bar(2, fcolor(black))
		
	gr export "AppendixFigure_4B.pdf", replace
	
	
	local i "att_oagv2d"
	g `i' =  `i'_18 -`i'_13
	replace `i' = 0 if `i' == .
	
	** Panel C. Attacks: Suitability
	graph bar `i', over(suit_median) ytitle("") ///
		graphregion(fcolor(white))  bar(1, fcolor(black)) bar(2, fcolor(black))
		
	gr export "AppendixFigure_4C.pdf", replace
	
	** Panel D. Attacks: Probability of substitution
	graph bar `i', over(pr_PNIS_50) ytitle("") ///
		graphregion(fcolor(white))  bar(1, fcolor(black)) bar(2, fcolor(black))
	
	gr export "AppendixFigure_4D.pdf", replace	
		
		
	
************
// Figure A.5. Evolution of coca prices in Colombia
************
	use "${Data}Data_PricesCoca.dta", clear

	two (connect Leafs year , yaxis(1) msize(small) msymbol(square) mcolor(black) lcolor(forest_green)) ///	
		(connect Paste year , yaxis(2) msize(small) msymbol(square) mcolor(black) lcolor(forest_green)) ///
		(connect Hydrochloride year , yaxis(2) msize(small) msymbol(square) mcolor(black) lcolor(forest_green)), ///
		legend(label(1 "Leafs") label(2 "Paste") label(3 "Hydrochloride")) ///
		xlabel(2005(3)2018) ///
		ylabel(1500(500)3000) ///
		ytitle("Price per kg") graphregion(fcolor(white)  /// 
		lcolor(white) ifcolor(white) ilcolor(white)) plotregion(fcolor(white) lcolor(white) ifcolor(white) ilcolor(white))
		
	gr export "AppendixFigure_5.pdf", replace	
	
************
// Figure A.6. Differential effects in protected areas
************
	* Load dataset
	use "${Data}MainData.dta", clear

	rename coca_protect_area coca_protect_area1 
	rename coca_No_protect_area coca_protect_area2
	keep coca_protect_area1 coca_protect_area2 coddepto year muni_code coddepto pr_PNIS_announce suit_announce iddeptoy announcement suitability pr_PNIS
		
	reshape long coca_protect_area, i(muni_code year) j(Type)
		
	g Protected = Type == 1
	g pr_PNIS_announce_Protected = pr_PNIS_announce * Protected
	g suit_announce_Protected = suit_announce * Protected
		
	g announce_Protected = announcement * Protected
		
	egen id = group(muni_code Type)
	
	** Panel A. Suitability
		* Plot
		forval v = 1/8 {
			g Prot_suitability_1`v' = Protected * suitability * (year == 201`v')
			g Prot_prPNIS_1`v' = Protected * pr_PNIS * (year == 201`v')
		}
		
		reghdfe coca_protect_area Prot_suitability_11 Prot_suitability_12 Prot_suitability_14 Prot_suitability_15 Prot_suitability_16 Prot_suitability_17 Prot_suitability_18 , absorb(id i.Protected#i.coddepto#i.year c.suitability#i.year) vce(cluster muni_code)
		preserve
			estimate store coef 
			parmest, norestore
			rename t tstat
			g year = 2010 + _n
			replace year = year + 1 if _n > 2
			set obs 8
			replace year = 2013 if _n == _N
			replace estimate = 0 if year == 2013
			replace min95 = 0 if year == 2013
			replace max95 = 0 if year == 2013
			replace parm="`parm_v'" if parm==""
			sort year
			
			twoway (connected estimate year , lc(gs0) ) ///
			(rcap min95 max95 year , lc(gs12) fc(gs12)), ///
			xlabel(2011(1)2018) ///
			xtitle("Years") ///
			ytitle("Coefficient") ///
			yline(0, lc(red) ) ///
			ylabel(-0.8(0.2)0, grid ) ///
			graphregion(color(white)) /// 
			legend(off)
			
			gr export "AppendixFigure_6A.pdf", replace
		restore		
		
	
	** Panel B. Probability of PNIS
		reghdfe coca_protect_area Prot_prPNIS_11 Prot_prPNIS_12 Prot_prPNIS_14 Prot_prPNIS_15 Prot_prPNIS_16 Prot_prPNIS_17 Prot_prPNIS_18 , absorb(id i.Protected#i.coddepto#i.year c.pr_PNIS#i.year) vce(cluster muni_code)
		preserve
			estimate store coef 
			parmest, norestore
			rename t tstat
			g year = 2010 + _n
			replace year = year + 1 if _n > 2
			set obs 8
			replace year = 2013 if _n == _N
			replace estimate = 0 if year == 2013
			replace min95 = 0 if year == 2013
			replace max95 = 0 if year == 2013
			replace parm="`parm_v'" if parm==""
			sort year
			
			twoway (connected estimate year , lc(gs0) ) ///
			(rcap min95 max95 year , lc(gs12) fc(gs12)), ///
			xlabel(2011(1)2018) ///
			xtitle("Years") ///
			ytitle("Coefficient") ///
			yline(0, lc(red) ) ///
			ylabel(-40(10)10, grid ) ///
			graphregion(color(white)) /// 
			legend(off)
			
			gr export "AppendixFigure_6B.pdf", replace
		restore		

		
